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The mixed spin 3- spin 3/2 Ising model has been simulated using standart 
CG ' 

. , and cooling algorithms on cellular automaton (CA) . The simulations have been 

-t— > 

made in the interval — 6 < D < 6 for J = 1 for the square lattices with periodic 
£^ . boundary conditions, including an anisotropy term D on spin-3. The model has 

o 

o 

compensation points through D / J = 2 line. The values of the critical exponents 

^ . ( p, a, /3 and 7) are estimated within the framework of the finite-size scaling 

O 

f"-> theory for selected D values (—2, 0, 1, 2 and 4). The estimated values are 

compatible with the universal two dimensional Ising values. 

CO 

I . INTRODUCTION 

The mixed spin Ising model has been preferred as a simple model to study 

ferrimagnetism which plays the fundamental role in the area of molecular based 

magnetic materials. In the ferrimagnetic material, the existence of compensation 

point makes some technological applications possible such as the thermomag- 



netic recording [1,2]. The various mixtures of spins have been studied with the 
simulation and the numerical methods as effective field theory [3 — 5], Green's 
function method [6] , mean- field approximation [7 — 9] , Bucndia and Novotnys 
conjecture [10], density functional theory [11] and Monte Carlo [12 — 20]. How- 
ever the mixed spin 3 - spin 3/2 Ising model has not been studied with any 
method yet. At the recent paper concerning the cyano-bridged terbium(/_ZT) 
- chrornium(77T) heterobimctallic assembly also attracted attention to absence 
of the study about spin 3- spin 3/2 mixture [21]. Therefore, the mixed spin 3- 
spin 3/2 Ising model is investigated via Cellular Automaton in this paper. 

In the previous papers, the Creutz cellular automaton (CCA) algorithm and 
its improved versions have been used successfully to study the properties of 
the critical behaviors of the Ising model Hamiltonians [22 — 32] . The CCA 
algorithm, which was first introduced by Creutz [33] , is a microcanonical al- 
gorithm interpolating between the conventional Monte Carlo and the molecular 
dynamics techniques. Our previous studies showed that the cooling algorithm 
improved from the Creutz Cellular Automaton algorithm are effective to study 
the phase space and the critical behavior of the Ising model [34 — 39] . 

The aim of this study is to examine the critical behavior of the mixed spin 
3- spin 3/2 Ising model. For this purpose, the order parameters (M, m a , mb, 



q a , qb), the susceptibility (%), the internal energy (Hj) and the specific heat (C) 
were calculated on square lattice L x L of linear dimension L = 32, 64, 128, 160 
and 256 with periodic boundary conditions. The {kTc / J , D) phase diagram in 
the interval — 6 < D < 6 for J = 1 was obtained using the temperature value 
determined from the maxima of the susceptibility (x) and the critical exponents 
were estimated by analyzing data within the framework of the finite-size scaling 
theory. The outline of this paper is as follows: In Section 2, the model and the 
formalism are given. In Section 3 the results and discussions are presented. In 
the final Section, a brief summary will be given. 

II . Model 

Two interpenetrating square sublattices occupied by the inequal spins com- 
posing the mixed-spin model. The Hamiltonian is given by 

H I = jJ2^S 1 +Dj2sf (1) 

where the sums are over all the nearest neighbors. J represents the bilinear 
interaction between spins <7i = ±3/2, ±1/2 on sites of sublattice A and Si = ±3, 
±2, ±1, on sites of sublattice B, D is a single ion anisotropy term on S = 3. 
The total energy H, which is conserved for the microcanonical algorithm, is 
given by 



H = H I + H K (2) 

where Hi is Ising energy defined in equation (1) and Hk is kinetic energy. The 
kinetic energy Hk is an integer, equal to the change in the Ising spin energy 
for any spin flip. For a site belonging to sublatticc A to be updated, its spin is 
changed to one of the other five states with 1/5 probability. If a site belonging 
to sublattice B to be updated, its spin is changed to one of the other three 
states with 1/3 probability. If this energy change is transferable to or from the 
kinetic energy variable of the site, considering the total energy H is conserved, 
then this change is done and kinetic energy is appropriately changed. Otherwise 
the spin is not change [34 — 37] . The simulations are carried out on the square 
lattice with L = 32, 64, 128, 160 and 256 for periodic boundary conditions. 

At the cooling algorithm, the simulation consists of two parts, the initializa- 
tion procedure and the computation of the thermodynamic quantities. In the 
initialization procedure, firstly, the sublattice A and B are decorated with S 
and <7, considering antifcrromagnetic interaction (J > 0). The initial configura- 
tion for the cooling algorithm is obtained adding energy to the kinetic energy 
of 40% of the spins on the every sublattice for getting the disordered phase (P) 
at high temperature. The initial configuration is run during the 20.000 Cellular 



Automaton time steps. During the cooling cycle, the kinetic energy (Hfc) of 
2% of the spins become equal to zero. Instead of resetting the starting con- 
figuration at each energy, it was used the final configuration at a given energy 
as the starting point for the next at cooling algorithm. The computed values 
of the thermodynamic quantities (the order parameter (M), sublattice order 
parameters (m a , TUb), the susceptibility (x, X a > Xb)> the Ising energy (Hj ) ,the 
specific heat (C )) and the Binder cumulant (Ul) are averages over the lattice 
and over the number of time steps (2.000.000 ) with discard of the first 100.000 
time steps during the cellular automaton develops [35]. 
Thermodynamic quantities are calculated by 



ma= N ^ S h m b = jj J2 a i ( 3 ) 

i(even) i(odd) 



I m„ — 771b I , x 

M = l - -± (4) 



_ (AP)-(M) 2 N (ml) - (m a ) 2 N (m 2 b ) - (m b ) 2 

X kT , X a 2 kT , Xb 2 kT 10J 



U I = (jJ2Si(T j + Dj2s 2 )/U (6) 



where Uq is the ground state energy at kT/J = 0. 

The phases are determined as 

Paramagnetic (P): m a = mb = 

Ferrimagnetic (/): m a =/= mb ^ 

III. Results and Discussions 

The ground state phase diagram in (D, J) is illustrated in figure 1 at the 
absolute zero temperature (kT/J = 0). It is seen that mixed spin 3- spin 
3/2 Ising model has ferrimagnetic (/) and paramagnetic phase (P). Through 
J = 1, the mixed spin system has ferrimagnetic /i phase as the sublattice A and 
B occupied by 5 — 3 and a = —3/2, respectively, in the interval —6 < D < 1.2. 
In the intervals 1.2 < D < 2 and 2 < D < 5, there are ji phase f^ ferrimagnetic 
ground state phases, respectively. 5 = 2 and a = —3/2 are located on the 
sublattice A and sublattice B creating ferrimagnetic ji while sublattices are 
decorated with 5 = 1 and a = —3/2 in f^. Both fy and P phases appears in 
the interval 5 < D < 6. For D < 6, there is paramagnetic phase (P). Through 
D/J = 1.2 phase separation line, 5 = 3 or 2 occupy the sublattice A while 
the sublattice B is decorated by a = —3/2. Similarly, sublattice A is occupied 



evenly by S — 2 or 1 while sublattice B is decorated by a = —3/2 through 
D/J = 2 phase separation line. So the generated values of S are limited to 
smaller values with increasing D. 

In figure 2, the temperature dependence of the sublattice order parameters 
are illustrated for selected D = I, 2 and 4 values representing ferrimagnetic 
regions for J = 1. As it is seen in figure 2(a), the mixed spin system exhibits 
second order phase transition (/ <— P). Both sublattice order parameters (m a 
and nib) appear continuous while the sublattice susceptibilities (x a and \ b ) 
have a characteristic peak at the same critical temperature for each D value. 
At D = 2, the sublattice order parameter m a is equal to 1.5 while the sublattice 
order parameter nib is equal to —1.5 cancelling each other for < kT / J, and 
hence there is a compensation point as Bobak et.al defined in the mixed spin-1 
and spin-3/2 ferrimagnetic system [40]. The sublattice magnetizations (m 
and mb) are illustrated for selected J = 1, 2 and 5 values through D/J — 2 line 
(figure 3). By this way it is verified that the sublattice order parameters ( m a 
and nib) cancel each other for D values through the D/J = 2 line. 

The lattice order parameters (M) and Ising energy (Hi) also appear con- 
tinuous (figure 4(a) and (b) ). The susceptibility (\) and the specific heat (C) 
have a characteristic peak which shifts to higher temperatures with increasing D 



value (figure 5(a) and (b)). It is obviously seen that single- ion anisotropy term 
D doesn't influence the phase transition type which is seen as of the second 
order. The simulations show that the mixed spin 3- spin 3/2 Ising system has 
second order phase transition from paramagnetic (P) phase to ferrimagnetic (/) 
phase (/ <r- P) using cooling algorithm in the interval —6 < D < 6. For D > 6 
, Ising energy doesn't exhibit any sign regarding a phase transition. 

In figure 6, the phase diagram is shown on the (kTc/J, D) plane through 
J = 1 line estimating the critical temperatures from the susceptibility maxima. 
In the interval — 6 < D < 5, the standart and the cooling algorithms give P — > 
/ second order phase transition at the same critical temperatures. However, 
standart algorithm still gives /->P second order phase transition for 5 < D < 
6, while there is no phase transition according to cooling algorithm results. The 
reason of this difference is occurrence of f$ and P phases as ground states in the 
interval 5 < D < 6. Thus, standart algorithm results initialized form the order 
state can exhibit /3 ground state phase while cooling algorithm results show P 
phase as ground state. 

The critical behavior of the quantities has been investigated within the 
framework of the finite size scaling theory for selected D values (—2, 0, 1, 2 
and 4) through J = 1 line. In here, the figures are represented for D — 1 which 



has /i ground state. The infinite lattice critical temperature Tc (oo) is obtained 
from the intersection point of the Binder cumulant curves for the different lat- 
tice sizes. As it is seen in figure 7(a), the Binder cumulant curves intersect at 
Tc(oo) — 4.61. The critical exponent v can be estimated using the finite-size 
scaling relation for the Binder cumulant given by 

U L = 1- -^ 4r 8 

3(M 2 } 2 

where reduced temperature is e = (T — Tc(oo))/Tc(oo) . The scaling data 
for the Binder cumulant lie on a single curve for universal correlation critical 
exponent value v = 1 of two dimensional Ising model (figure 7(b)). X'c(oo) is 
also obtained by the finite-size scaling theory using following relation 

T c (oo) = T c (L) + aL- 1 ^ (9) 

Extrapolation of the susceptibility peak temperatures to 1/ ' L x l v — > gives us 
Tc(oo) value. Both from the intersection temperature value of the Binder cumu- 
lant curves and the extrapolation of the lattice susceptibility peak temperatures 
for different lattice sizes give Tc(oo) = 4.61 at D — 1 (figure 7(c)). 

In order to investigate the critical behaviors of the lattice order parameter 
(M), susceptibility (\) and the specific heat (C), the following finite-size scaling 



relations are used near the lattice critical temperature. 

M = L~^ v X{eL 1 ' v ) (10) 

kT x = IT'l v Y(eL 1 f v ) (11) 

C = L a / v Z{eL 1 / v ) (12) 

The infinite lattice critical behaviors must be asymptotically reproduced for 
large x — eL 1 ^, that is, 

X(x) = Bx fi (13) 



Y(x) = Cx-~* (14) 



Z{x) = Ax- a (15) 

The finite-size scaling plots of the order parameter are exhibited in figure 
8(a) at D = 1. For f3 = 0.125 and v = 1, data lie on a single curve both below 
Tc(oo) and above Tc(oo). For T < Tc(oo), the data of order parameter M are 
in good agreement with the universal value of ft — 0.125. For T > Tc(oo), the 

10 



data fit to the straight line with slope = j3 = 0.875 according to the equation 
(14). In figure 8(b), the finite size scaling plots of the data for the susceptibility 
are illustrated together the straight lines describing the theoretically predicted 
behavior for large x (equation (15)). The scaling of the susceptibility data agrees 
with the asymptotic form with 7 = 7' = 1.75 and v = 1 for both sides of the 
Tc (00). The specific heat of an infinite lattice is well described by 

C/k = Ae~ 0: + b ± (16) 

where b^ expresses the nonsingular part of the specific heat [41]. The finite-size 
scaling plots of the singular portion of the specific heat (C/k — b ) arc shown 
in figure 8(c) and figure 8(d) for the selected D = 1 value. The data for the 
specific heat scale well with the values of a = 0.03 at T < Tc(oo) and a ' = 0.01 
at T > T c (oo). 

At D — 2 where m a and m^ compensate, the infinite lattice critical temper- 
ature Tc(oo) is obtained as 2.96 ±0.01 from the intersection point of the Binder 
cumulant curves (figure 9(a)) and from the extrapolation of the susceptibility 
peak temperatures to l/L 1 /" — > (figure 9(c)) for the different lattice sizes. On 
the other hand the universal correlation critical exponent value is obtained as 
v = 1 from the scaling of the Binder cumulant (figure 9(b)). In figure 10, the 
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critical exponents are estimated for D = 2 where the m a and mj cancel each 
other. They are estimated as v = 1, ft = 0.125, ft' = 0.875, 7 = 7' = 1.75, 
a = 0.05 and a ' = 0.04 within the framework of the finite-size scaling theory 



are compatible with theoretical values. 

The estimated values of the critical exponents (a, ai, ft, ft/, 7, 7/ and v), 

the infinite critical temperatures as well as b^ are given in table 1 for selected 

D values (—2, 0, 1, 2 and 4) through J = 1 line. As it is seen from table 1 that, 

the estimated critical exponents are in good agreement with theoretical ones 

(y = 1, a = a' = 0, ft = 0.125, ft' = 0.875 and 7 = 7' = 1.75 ). 

Table l.The estimated values of the the critical exponents 

T^(oo) 

T*(oo) 
b- 

b+ 
v 
a 

a' 
/3 
/?' 

7 

i 

III. SUMMARY 



D= 2 


D = 


D = l 


D = 2 


D = 4 


7.46 ±0.01 


5.99 ±0.01 


4.61 ±0.01 


2.96 ±0.01 


1.46 ±0.01 


7.46 ±0.01 


5.99 ±0.01 


4.61 ±0.01 


2.97 ±0.01 


1.46 ±0.01 


0.01 


0.06 


0.01 


0.015 


0.05 


0.02 


-0.07 


-0.03 


0.050 


0.00 


1 


1 


1 


1 


1 


0.010 


0.008 


0.030 


0.050 


0.030 


0.010 


0.020 


0.010 


0.040 


0.050 


0.125 


0.125 


0.125 


0.125 


0.125 


0.875 


0.875 


0.875 


0.875 


0.875 


1.75 


1.75 


1.75 


1.75 


1.75 


1.75 


1.75 


1.75 


1.75 


1.75 



The two dimensional spin 3- spin 3/2 mixed spin Ising model has been sim- 
ulated using cooling algorithm on the cellular automaton (CA) for the square 
lattices in the interval — 6 < D < 6. The simulations show that the spin system 
with antiferromagnetic interaction has ferrimagnetic (/1 , f2, J3) and param- 
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agnctic (P) phases at absolute zero temperature (kT/J = 0). The existence 
of single-ion anisotropy term D does not effect the phase transition type but 
the phase transition temperatures shift to the low temperatures with its in- 
creasing value. As it is recognized from the {D, J) phase diagram that the 
generated values of S are limited to smaller values with increasing values of 
D. An other result is occurrence of the compensation points through D/J = 2 
line at kT/J = 0. For selected D values (-2, 0, 1, 2 and 4) through J = 1 
line, the estimated values of the critical exponents which are independent from 
D value are in a good agreement with universal values of the two dimensional 
Ising model. The standart and the cooling algorithm are suitable procedures 
for simulating mixed spin systems. 
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Figure Captions 

Figure 1. The ground state phase diagram of the mixed Spin 3- -spin 3/2 
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fcrrimagnetic Ising model. / and P indicate ferrimagnctic and paramagnetic 
phases, respectively. On the phase separation line between ferrimagnetic /i and 
/2 phases there is fcrrimagnetic f^ phase. 

Figure 2. The temperature dependence of the (a) sublattice magnetizations 
(m a , mi,) and (b) the susceptibilities (x a , Xb)- 

Figure 3. The temperature dependence of the sublattice order parameters 
(m a and mt) through D/J = 2 line for J = 1, 2 and 5. 

Figure 4. Temperature dependences of (a) Lattice order parameter (M) and 
(b) Ising energy (Hi) . 

Figure 5. Temperature dependences of (a) lattice susceptibility (x) and (b) 
The specific heat (C/k) . 

Figure 6. Phase diagram of the mixed spin 3- spin 3/2 Ising model on 
(kTc/J, D) plane for J = 1. Thick line and square symbol indicate the standart 
and the cooling algorithm results, respectively. Vertical thin lines are used for 
separating ground state regions. 

Figure 7. For D = 1 and J = 1; (a) Temperature estimations from the 
intersection point of the Binder curves, (b) The estimation for the critical ex- 
ponent v from Binder cumulant scaling and (c) Temperature estimations from 
the linear extrapolation. 
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Figure 8. For D = 1 and J = 1; finite-size scaling plots of (a) the lattice 
order parameter (M), (b) the susceptibility (x), (c) the specific heat (C) at 
T <T C and (d) the specific heat (C) T >T C . 

Figure 9. For D = 2 and J = f; (a) Temperature estimations from the 
intersection point of the Binder curves, (b) The estimation for the critical ex- 
ponent v from Binder cumulant scaling and (c) Temperature estimations from 
the linear extrapolation. 

Figure 10. For D = 2 and J = 1; finite-size scaling plots of (a) the lattice 
order parameter (M), (b) the susceptibility (x), (c) the specific heat (C) at 
T < T c and (d) the specific heat (C) T>T C . 
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